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Abstract 

Two finite difference schemes for the numerical treatment of the von Neumann equation for the (2+l)D 
Dirac Hamiltonian are presented. Both utilize a single-cone staggered space-time grid which ensures a 
single-cone energy dispersion to formulate a numerical treatment of the mixed-state dynamics within the 
von Neumann equation. The first scheme executes the time-derivative according to the product rule for ” bra” 
and ”ket” indices of the density operator. It therefore directly inherits all the favorable properties of the 
difference scheme for the pure-state Dirac equation and conserves positivity. The second scheme proposed 
here performs the time-derivative in one sweep. This direct scheme is investigated regarding stability and 
convergence. Both schemes are tested numerically for elementary simulations using parameters which pertain 
to topological insulator surface states. Application of the schemes to a Dirac Lindblad equation and real- 
space-time Green’s function formulations are discussed. 

Keywords: Dirac equation, vonNeumann equation, finite-difference scheme, staggered grid, fermion 
doubling, topological insulator 


1. Introduction - preliminaries and definitions 

1.1. The Dirac equation and numerical schemes 

In context with graphene and, more recently, topological insulator surface states, the Dirac equation 
has received renewed interest within the physics community. Introduced by P. A. M. Dirac in 1928, the 
(3+l)D Dirac equation traditionally has been a model in high-energy physics to describe relativistic spin-1/2 
particles and, as a field theory, is a key ingredient to the standard model of elementary particle physics. m 
In condensed matter and atomic physics its importance lies in the non-relativistc limit providing the spin- 
orbit interaction, which is instrumental to an understanding of atomic spectra and represents the foundation 
for the field of spintronics.® In electronic structure calculations the full Dirac equation has been used to 
describe inner shell electrons. 0® While the (1+1)D and (2+l)D Dirac equation, allowing a two-dimensional 
representation of the Clifford algebra, in the early days of quantum physics was used for simplicity in formal 
model studies, condensed matter physics has recently identified physical realizations of the (2+l)D Dirac 
equation in form of low energy electronic excitations in graphene and topological insulators (TI). [9HT4] 

With the interest in graphene and TIs as components in electronic and spintronic devices efficient schemes 
for the numerical solution of the (2+l)D Dirac equation have become desirable. Numerical approaches for 
solving the Dirac equation have taken several approaches. For investigations of relativistic electrons in 
atomic physics, a (2+l)D FFT-split-operator code was usedlTS:, TBJ. In such an approach, fast Fourier 
transformation between the space and momentum representation is used to compute the Dirac propagator. 
The computational effort scales like 0(N\nN) where N is the number of grid-points. An efficient code 
using operator splitting in real space was introduced for the (3+l)D case. jT71 It leads to the highly efficient 
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operations count of O(N). Alternatively, approaches using finite-difference approximations to the first- 
order space-time derivatives have been used and have foremostly been applied in lattice quantum chromo 
dynamics (QCDl.[l8] Traditionally they were plagued by the fermion doubling problem associated with 
centered difference approximations to first-order space- and time-derivatives. [H5] We note that such a 2d 
finite-difference scheme for the Dirac-Poisson system was presented recently. p?0] 

Recently we have identified a class of staggered grids on which fermion doubling can be avoided in 
arbitrary space dimensions, with explicit schemes given for the (1+1 )D, (2+1 )D, and (3+l)D case. ITl. 
This class of schemes scales linearly in the number of grid points and allows a gauge-invariant formulation of 
electromagnetic fields on the lattice. Moreover, this grid can be used in any space-time formulation involving 
the Dirac operator, such as the density matrix and Green’s function approach, while yielding a single Dirac 
cone. 

In this work we develop two numerical schemes for the mixed-state time evolution under the Dirac 
Hamiltonian H in (2+l)D , which we term the Dirac-von-Neumann equation 


ihp = Hp — pH = [H, p] . (1) 

It describes the coherent mixed state dynamics of Dirac fermions which is useful for simulations of quantum 
transport on TI surfaces near the Dirac point. Both schemes to be presented in the following sections utilize 
the staggered grid of Ref. [22] to define centered differences over a single lattice spacing, thereby eliminating 
the very source for fermion doubling from the start. The first scheme (”bra-ket”) scheme, introduced in Sect. 
[2j treats the time derivative of the density operator p = {4>i\ within the product rule as indicated 

by the commutator on the r.h.s. of Eq. [I] apply H from the left and therewith propagate the ket in time, 
apply H from the right and therewith propagate the bra in time, and then form the difference. The second 
(’’direct”) scheme treats the time-derivative of p in one step and is introduced in Sect. [3l with further details 
given in the Appendix. Numerical examples, one for each scheme, are given in Sect. [4] Finally, we briefly 
discuss the extension of this approach to the Lindblad equation and space-time formulations of the Green’s 
function method for the Dirac Hamiltonian in Sect. [5] and close with Summary and Outlook. 

1.2. The continuum formulation of the von Neumann equation for the (2+l)D Dirac Hamiltonian 

We consider the effective model Hamiltonian which accounts for the energy spectrum of TI surface states 
near the Dirac point [lOQ 

H = v(axP) \ z +m(X,Y,t)a z + V{X,Y,t)l 2 (2) 

Oi denote the Pauli matrices. Using the abbreviatiorj^] 

d± = v ( P v ± iP x ) 


and (omitting space-time arguments) 

m± = V ± to 

this Hamilton operator takes the form of a 2 x 2 matrix 

( m+ d+ \ 

H ~\d. m_ ) ■ 

Inserted into the von Neumann equation ([!]) one obtains a set of first-order partial differential equation in 
space and time for the density operator elements conveniently written as the 2x2 matrix identity 

ifr— ( Pn pl2 \ = ( m + Pl1 “ Pu m + + d +P2i - Pi2d- m + pi2 - Pi2m_ + d+p 2 2 - Pi\d+ \ 

dt\P2i P22 ) \ m-p2i-p2im + +d-p 11 -p 2 2d- m_p 2 2 - p22m- + d-p 12 - p2id+ ) 


1 Note that the TI representation differs from the standard representation used in Ref. m, however, shifting P x —>■ P y and 
P y —> —P x converts the latter to the former. 

2 For now, the electromagnetic vector potential is set equal to zero. 
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Note that, as in the original von Neumann equation, Hamilton matrix operator elements to the right of the 
density operator matrix elements act to the left (and vice versa). The two-component nature of the spin-1/2 
Dirac fermion suggests this 2x2 form. Choosing a continuous space representation for the orbital degrees 
of freedom one arrives at 

ihp(r,r') = (3) 

f (m+( r) - ra + (r'))pn(r, r') + d + p 2 i(r, r') + d'_p 12 {r, r') (m+(r) - m_ (r'))pi 2 (r, r') + d + p 2 2 (r, r') + d' + p n { r, r') 
V (w-(r) - m + (Y'))p 2 i(v,r') + d-pu( r,r') + d'_p 2 2 (r,r') (m_(r) - ra_ {r'))p 22 (r, r') + <9_pi 2 (r,r') + d' + p 21 (r,r') 

(4) 

Here we omit the time variable t, common to all terms, for brevity, and use the real-space versions of the 
abbreviations defined above, be., 


d± = — (d y ± id x ) , &± = — ( d y > ± id x >) 

and (including arguments) 

m±(x,y,t) = V(x,y,t) ± m{x,y,t ) . 

Note that in the real-space representation we have 

{$3 I ( P y ± iP x ) I x,y) = (( P y =F iP x ) ^j | x,y) = (x,y \ {P y =F iP x ) %)* = -d± i & j (x,y)* . 


1.3. Placement onto a space-time grid 

The task at hand now is to develop a space-time finite-difference approximation to Eq. [4] which avoids 
fermion doubling. This is achieved by using the staggered grid introduced in an earlier paper to accommodate 
the real-space density matrix in which upper and lower spinor component(s) are placed on two adjacent time- 
sheets. [23 The proper implementation is found by inspection of the density operator for a pure state spinor. 

ipi 


In (2+l)D one has a pure-state ket | US') = 


i>2 


and a ket-bra projector for the density operator 


p= l*X*l=(v£)(^ M) 


ipripl ip \ 

1p2lpl 1p2lp 2 J 


( 5 ) 


This shows that the Pauli indices of p , 1 and 2 respectively, take the face-centered rectangular ^>i-grid (Q i) 
and ip 2 -grid (Q 2 ) position. For a single-cone representation the latter are[22) 


V’i(j) with j e Q x = {(j x , jy,j 0 ~ 1/2), {j x + 1/2, j y + 1/2, j a - 1/2) | j v e = x,y,o} , 

V> 2 (j) with j e Q 2 = {{j x + 1/2, j y ,j 0 ), Ux,jy + l/2,jo)\jv€'Z-,v = x,y,o}. 

Here the time index is labeled o. Note that, for given time, ipi and ip 2 are placed onto adjacent time 
sheets, respectively, jo — 1/2 and jo- Each time sheet contains a rectangular face-centered spatial grid, 
symmetrically staggered relative to the ones on the two adjacent time sheets, such that symmetric difference 
quotients replace the respective partial derivatives on the grid. The grid spacings are denoted by A x . A y , 
and A 0 = vA t . We also introduce the associated grids 

Qi = {{jx, jy,3o), {jx + 1/2 ,jy + 1/2, jo) \j„ez,v = x,y,o} , 

02 = {{jx + 1/2, jy, jo + 1/2), {jx, j y + 1/2, jo + 1/2) I j v e T,V = x,y,o} . 

These grids are obtained from their non-barred counterparts by a forward time-shift by A t /2 and thus share 
the spatial positions with the former. These will be the grids on which we place mass, scalar potential, and 
all partial derivatives. Common to all four grids is their symmetry in space-time, as required in a covariant 
theory. Since most dynamic problems are formulated as an initial-value problem in time, however, it is useful 
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to view each of the two grids as a set of time sheets. Each time sheet in turn consists of a rectangular-, for 
A x ^ A y , or cubic-face-centered spatial grid, for A = A x = A y , which in units (A x ,A y ) is characterized by 
the set of relative positions ( j x ,j y ) £ Z 2 U(Z + |) 2 . 

The matrix elements pij (and the vector potential to be introduced later) live on the grids Qi 

Pij(r,t-,r',t') -> PijiGiiGj) , i,j = 1,2. 


Note that t = t' leads to the need for two adjacent time sheets (e.g., j 0 — 1/2 and j 0 ) for the representation 
of all matrix elements of p{t) on the lattice: while pn and pii each are placed on a single time sheet, p\i 
and P 21 require the use of two adjacent time sheets. This leads one to the definition of the trace of p using 
two adjacent time sheets of the grid: 

Definition 1. The trace of p at a given time (j 0 — 1/2 ,j a ) on the grid, 


Tr {p} 


Tr{p 11 (^ io - 1/2) ;^ 0 “ 1/2) ) + P22(# ) ;# ) )} 


E 

(jx ,jy ) 


Pll ( jxijy) jo 


1 . . . 
2 1 JxiJm Jo 


2) + P22 (jx + 


1 

^ > Jy)JoiJx 


) JytJo) 


( 6 ) 


is defined as the sum of the trace of p\\ with respect to the lattice sites Q^° and the trace of P 22 with 


,Uo) 


The summation runs over all lattice sites of the respective time sheet, 


respect to the lattice sites Q^ 

(jx,j y ) £ Z 2 U(Z + \) 2 . 

We extend this definition to (j x ,j y ) lattice sums for non-diagonal density matrix elements 
Tir{pij{Q;Q )} = \Pij{jx,jy>jt'ij x ijy>jt)\ \j! c =jx+Ajx,jy=jy+Ajy > 

{Jx j jy) 

with constant relative indices Aj x and Aj y . 

Definition 2. Centered spatial difference operators are defined on the barred grids as follows 
D x f{Qi) ■■ D x f \jl Jy = ~ Ux,jy,jt ) € Qi , 


(7) 


Dyf{Qi) : D v f | j a .,j y — » (fjt,jy + 1/2 fjx,jy- 1/2)1 (, jxijyijt ) £ Qi 


D±f(Qi) : D±f 


\3t 


A a 

= Dyf 


±iD x f 


\3t 


and 


1 


dt+1/2 


D om): D of ^ tj =-(f^ y 


_ ft- 1/2 
J 3x iJy 


1? (jxi 3 y> jt) £ Qi ? 


)? (jxi jyi jt) £ Qi 


Furthermore, we will use the time-average operation 

\3t 


r/(ft) : Tf |>; A = U,M,) € ft ■ 


’ 3 x i 3 y 

Formally these operators are defined on the barred grids Qi, however, the execution involves density 
matrix elements on the unbarred grids. Note that when applied to the density matrix, these operators may 
act on the first (bra) or second set of indices (ket). 

Partial derivatives, scalar potential and mass terms, respectively, are implementing into Eq. ([ 3 ]) by the 
substitutions 

Q vh 
d± —> —D± , 


m±(x,y,t) [V(Q i )±m(Q i )\T=2ivhM±(Q i )T, M± = - 


rn± 


2ihv 

Within a sum over all lattice sites ( j x ,j y ) £ Z 2 U (Z + \) 2 a spatial difference operation {D' k , k=x,y) 
performed on the second set of indices i' of the density matrix p{i\ i') is equivalent to minus the derivative 
(— D k ) performed on the first set i. 
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Corollary 1. Integration by parts on the grid: Let ( jx,jy ) € Z 2 U (Z + l) 2 . Let D k and D' k , k = x,y as in 
Def. i respectively, denote spatial difference operators with respect to the first and second set of indices of 
Pij. Then, under periodic or zero spatial boundary conditions for p and with Def. [7] 

Tr {D' kPl fig-Q')} = Y D' kPij (Q-g') = - Y D kPii & Q') = - Tv{D kPij (g-Q')} , fork e {*, y} . (8) 

(jx ijy) (jxijy) 

Proof. We give the proof for i = 1, j = 2, k = y 


(jx ijy ) 


E 

(jx ij y ) 

1 


P 19 ... i P 12 .... ] 


Ayl 


Y p^r hjo 


1 E^r l;j ° 


A v i ^ rL 2 j*,j y ;jx'J y ' + h Ayl ^ i 

" Cj*,j») " (j*,j„) 

1 

A/ 


E 


(j*->j*+i/2„j y ->j y +i/2) 


Pl2 jx-l/2,jy-l/2-,j x ,-l/2,j 


1 

A 


E jo- 2 ;jo 

W2 j<r-l/2,j»+l/2;j a! /-l/2,j v , 

^ (j* —>j* + l/2,jy—>j y —1/2) 


- E 

(jx ijy ) 


2 ;SrV 


(9) 


We have used that A y / = A y and that, under the trace (j x / — j x ) and (j y / — j y ) are constant. Note that a 
combined shift (carried out in the third step) of j x and j y is required to stay on the proper sub-grid f/,;. This 
mixes corner and face-center positions of a given time sheet. However, since the sum is to be taken over 
all grid points this ’’integration by parts” rule holds under the trace. Note also how the barred grids used 
for the definition of center position of the spatial derivative are related to the index of the density matrix 
elements. Other cases for i, j , and k can be shown in the same fashion. □ 


2. The bra-ket scheme 


The bra-ket time-propagation scheme is obtained directly by adapting the scheme for the pure-state Dirac 
equation, applied to the bra- and the ket-side of the density operator. This corresponds to the interpretation 


p = Y^ k [i**x*ki + i**x**il = 


( 10 ) 


Definition 3. Within the representation Eq. ([2]) the pure state time-evolution of Ref. from initial time 
” — ” : jo — 1 / 2 , jo to final time ” + ” : jo + 1 / 2 , jo + 1 may be written as 


with the coefficients 


i>t = <#i + 2 

1 P 2 = 1^2 + Wl , 


a = 


A a 


■ Mi 


±~ M + 


a* = 1/a , 


P = ~ — 




— M_ 


-D+ , 


-D- , 


( 11 ) 

( 12 ) 


5 


Wo + M + 








living on integer time sheet j Q of grid Q\, and 


£+ M - 

±~ M - 


7* = 1/7 , 


5 = —- 


1 




D_ , 5* = — 


■ M_' 


-D 


+ > 


living on half-integer time sheet j Q + 1/2 of grid Q 2 ■ We have used ( M ±)* = —M± which holds for real-valued 
mass and scalar potential V. 


Using the short-hand notation p\y = (x,y \ p(t ) | x',y'), the progression by one full time step A 0 is 
executed as follows (initial (—) and final (+) time is indicated by the respective superscript): 


Pi-i'- 


P 1 - 2 ' — 

P‘1-1 '- 


P2-2'- 


This set of operations is followed by 


Pi+i'- a P\-i'~ + &P2-V- 

Pi-v+ = a ' P\-v~ + P' P\- 2 '~ 

Pl+2'- = 0+1-2'- + PP2~2'~ 
Pl- 2 '+ = l' Pl- 2 >- + P' Pl-v + 

P 2 +V- = IP 2 -V- + Ppl+V~ 
P2-V+ = a ' P2~V~ + P' P2-2'~ 

Pl+2'- = 1P2-2'- +^Pl + 2'- 
Pl-2'+ ~ l' P2-2'- + P' P2-V + 


Pl+1'+ — a Pl~l'+ + PP2~1' + — a ' Pl+1'~ + P' Pl + 2’~ 

Pl + 2’+ = &Pl-2'+ + PP2~2’ + = l’ Pl+2'“ + P' Pl+V + 

P2+V+ = 1P2-V+ + ^Pl+1'+ = a ' P2+1'- + P' P2+2'- 

P2+2'+ = lPl-2'+ + ^Pl+2' + = l' P2+2'- + P' P2+1' + 

Executed in this order the scheme is explicit and follows exactly the single-cone time-propagation of Ref. 
)22j applied to both sets of indices of the density matrix, with unprimed and primed operators acting on, 
respectively, bra and ket. As such all properties of the scheme for the pure-state propagation, such as 
stability, convergence, and spectral properties, are carried over to the scheme for the density matrix. 


Definition 4. With the abbreviations defined above, the single time-step progression under the bra-ket 
scheme is defined as 


Pi+v+ — 01 a! Pi-y- + aft 1 Pi-y- + P a> P 2 - 1 '- + PP' p 2 - 2 '~i (13) 

Pl + 2’+ = 017 ' Pi-2'- + Pi' p2-2>- + P' Pl+V+ > (14) 

P2+l'+ = l a '* P2-1'~ + lP'*P2-2'~ + PPl+V+ , (15) 

PPl+2' + ■ 
(16) 

Def. [4] extends the single-cone pure-state time-propagation scheme in Def. [3] to the one for mixed-states. 
For pure states, the latter reduces to the former. 

The following lemmas and propositions pertaining to the bra-ket scheme are a direct consequence of the 
properties of the pure state scheme Ref. |22l and the observation that the density operator at initial time, 
without loss of generality, can be written as as sum of pure-state contributions of the form Eq. [5j 


P 2 + 2 '+ = l'* (7 + PP) P 2 - 2 - + Pay* p 1 - 2 >- +5'*p 2 +y+ = 7 ( 7 '* + P'*P'*) P 2 - 2 '- +7 P 


/ * 

a p 2 - 
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Lemma 2. The bra-ket scheme conserves positivity and Hermiticity of the density matrix. 

Proof. Both properties are a direct result of the bra-ket scheme: Let, at initial time the density matrix 
Po = J2kTk I ^k){^k I be positive definite, i.e., 


($ I Po I $> = I tfPo&j) O' I $} > 0 


for an arbitrary element | $) of the Hilbert space. Abbreviating the one-step time evolution of a pure state 
Eqs. © and ^ by = K^~, 


K = 


a /3 
1 5a (7 + S/d) 


the bra-ket scheme propagates 

<$ I Po I $) -»• <$ I K Po rf I $) = <$' I Po I $') > 0, I $') = Jft | $) 


since positivity holds for p Q . Hence under the bra-ket scheme positivity is preserved for arbitrary time step. 
Preservation of Hermiticity is seen from ( Kp 0 K 1) = Kp a KY □ 

Note that the standard norm is not conserved within the pure-state scheme since K is not unitary. [22] 
Hence, although the time-evolution superoperator can be cast in Kraus form with a single Kraus operator 
(living on two adjacent time-sheets), we do not have K = l. (251 IM] 

Lemma 3. Under periodic, respectively, zero boundary conditions, the functional 

E° r = TT{p u (g[ j °- 1/2) ;G[ j °- 1/2) ) + P 22(Gi jo) ;Gi jo) ) + A 0 K[D_p 12 (^ 0 “ 1/2) ;^ o) )]}, (17) 

with r = ( r x = A a /A x ,r y = A 0 /A y ), is conserved under the bra-ket time propagation, /d/.s[7~1 to\10\ 

Proof. The proof of this lemma as well as the following two propositions rest on the fact that the density 
matrix at initial time may be cast into diagonal form and then, under the scheme, is placed onto the grid 
according to 


p=J2 ^ i * (fc) x* (fc) 


7 k 


,(fc) Ak)* 

1p\ 'lf\ 

Ak) Ak)* 

V>2 Vl 


4 k W 

v4 fe) vT 


E 


7 ?= 


4 k) (Si)4 k] (Gi)* 4 k) (Gi)4 k) (G 2 )* 


4 k] (G2)4 k) (Gi) 


4 k) (52)4 k \G2Y 

(18) 


where | are normalized to one, 0 < 7*, < 1 real, and X) fe 7fc = 1. Since each term in this sum i 


is 


propagated individually within this scheme, all the stability properties of the pure-state dynamics apply. 
The conserved functional for the individual pure-state contribution, E^ k \ and the definition of the trace on 
the grid in Eq. as a sum over all lattice sites of the respective time-sheet render the conserved functional 
for the density operator JYkTkE/^ = E® . □ 


Proposition 4. Let r x = 


r v = 


Ao 


with r x + r y < 1 (e.g. using r x = r y < 1/2). Then, the bra-ket 


scheme Eqs \13\ to 10 is stable and satisfies the estimate 
Tr{pn(S[ 


Uo - 1/2) _GU.-i/2) )+P22 ^ 


? b) ;G?° 


')}< 


El_ 

1 -r x -r v 


(19) 


for all time. 


Proof. Using the decomposition Eq. (181 stability can be shown term by term in the sum. 


.□ 
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Comment: Conservation of E® on the grid corresponds to the conservation of the trace of p in the con¬ 
tinuum limit. At the same time the conservation of the former within this scheme implies non-conservation 
of the trace when the latter is defined on the grid according to Eq. § This non-conservation is of order 
A 0 /A and thus can be adjusted by this ratio. 

From the analysis of the pure-state scheme we know that Prop. [4] is too restrictive. In fact, the stability 


condition r 2 + r~ < 1 holds for constant mass and potential. For the special case r = r x = r y this less 


restrictive stability condition reads r < 1/ \[2. This stability condition also holds for an arbitrary space- and 
time-dependent mass and potential. 

Definition 5. We define the functional 


Mp{°i 5 ) = 


2 ^ 

jxj y 


2 {jx, jy, jx, jy) + 


P 11 


Pll 2 {jx, jy T 1 ; jxi jy) (■ + Pu 2 {jx + „ , jy + 0 i jx 


Jo -3 


2 ’“ 


+ 

+ 


s <! Pu 2 


l 


l 


{jx “b 2 lis “f 2 ’ J*’ jy^ ^ Pll ” {jx T 2 ’ jxi jy 1 ) 


.Jy + 2 )} 

(20 


%<Pii 2 


{jx 2.Jy (j* 9 > Jy + 2’Z.Jy + l) 


Tr(p J 


22 ) 


= E 


Jx jjy 


PyzUxijy Wljxijy o) + r, 


^ P22{jxijy + 0 >jx> jy 0 ) f + P 22 {jx + „ . jy , jx n,jy) 


+ ^ (^{jx + 2 ’Jy>J*’Jy 2) 1 ^| ^22O'® + 2’J'«>J'*’J'w 2) 


( 21 ) 


+ Pmijx 2’Jy’J^’Jy 2) 1 ^22(j* 2.Jy.J®.Jy 2^ 


Again, the summation runs over all lattice sites of time sheet Q\ l , i.e., j x ,jy £ Z 2 U (Z + with 


ji = jo - o J 2 = Jo Jo £ 


Proposition 5. Tei r = r x = r y = l/y/2, hold in (13 1 -(16). XTien this scheme is stable and for all time 
satisfies 

Tr(pfr 1/2 ) + Tr( P r 2 ) < 2£° , I? 0 = ^ . (22) 


Proof. This inequality has been shown to hold for pure states based on the definition of an averaged norm. [22] 
In the present representation of the Dirac Hamiltonian Eq. ([ 2 ]) one defines for upper and lower spinor 
component ipi and fa, respectively, 


jj° 2 

2 

E+I,iv , v 'i" E+ Uv-h-^ 

Vi 

2_^ 

(jx,jy) 

2V2 2V2 


(23) 


and 


V> 2 ° 




>Jy 


V>2 


2 


2^2 


2V2 


(24) 


- E 

Ux jy) 

with j x ,j y £ Z 2 U (Z + 1/2 ) 2 . For a pure state the averaged norm can be recast into the expression (21) for 


pu on subgrid 2 . Similarly, using the averaged norm for components fa one arrives at (22) for subgrid 
Q 2 ° ■ Using the pure-state decomposition Eq. (18) , adding the 7 j weighted contributions for individual pure 
states verifies inequality ( 22 ) for a general density matrix. □ 


3 For compactness of notation the diagonal time-index is placed onto p as a superscript. 



































3. Direct scheme 


An alternative, direct scheme using the same placement of density matrix elements on the lattice formu¬ 
lated is formulated here. It allows centered difference quotients over one lattice spacing for an approximation 
of all first-order partial derivatives thus eliminating the source for fermion doubling. Rather than follow¬ 
ing the separate bra-ket application of the pure-state scheme, however, it is based on a simultaneous time 
propagation of the bra and ket indices thereby reducing the number of operations required per time step. 

The compact global form oft his numerical scheme may be formulated as an initial value problem: at 
initial time the four density matrix elements are stored on the two adjacent time sheets, Gi° 1 ^ 2 ' > and G 2 \ 
associated, respectively, with time j a — 1/2 for grid Gi and time j 0 for grid Qo- Replacing partial derivatives 
in Eq. [3] by difference quotients, as detailed in Sect |1.3[ on obtains for this direct scheme: 

Definition 6. 


D 0 pn(Q[ jo) ;g[ jo> ) = - M + (QV°> ) TpMG 


(jo). 




iUo).nUoY 


1 ^1 


- D +P21 (g[ Jo) -,g[ jo - 1/2Y ) - D'_ Pl2 (g[ j °- 1/2) -,g [ j ° y ), 


(25) 


D +P2 2(G[ Jo) -,Gi io) ) - D / +Pl gg[ Jo+1/2) -gi jo+1/2) ) , 


(26) 


D 0 P2i(g ( 2 lo+1/2) ;g[ J ° y ) = 


T(A+1/2) 


M-(g 2 

D_ m (^ 0+1/2) ;# +1/2) ' 


) - M + (g[ joY )} Tp 2 i(Gi j ° +1/2) ;G[ joY 


) - D'_p 2 2(G ( 2 3o) ;G[ joY ) , 


(27) 


and 


D 0 P22(G 2 


(ic + l/2).^bo + l/2)'- ) 


= M_(^ +1/2) ) - M_(^ +1/2) ')] Tp22(G { 2 jo+1/2) -,G [ 2 Jo+1/2y ) 

- D_p 12 (g?° +1/2) -G { 2 io+lY ) - D' + p 21 (Gi jo+1) ;G ( 2 jo+1/2Y ) ■ ( 28 ) 


Eqs. (26) and (27) are equivalent. The explicit form of the scheme for specific lattice sites is given in 


the Appendix. 

The key, as for the pure state case, is a combined staggering in space and time which allows the calcu¬ 
lation of symmetric difference quotients for each of the first-order partial derivatives. Note that the first 
order-difference operators, as well as the time-averaging operator, which formally are defined on grids Gi 
place the density matrix elements on the grids Gi (rather than Gi), so that the latter must be stored as 
matrices pij(Gi', Gj ), i,j = 1,2. In particular, spatial derivatives are performed on the time-sheet which is 
sandwiched between the initial and final time-sheet of the associated row or column which is propagated: 
spatial derivatives map g 2 ° +1 ^ 2 ' > onto Gi° +1 ^ 2Y , likewise is mapped onto G 2 "^- The numerical proce¬ 
dure follows a local time-sheet by time sheet propagation which allows full parallelization in its execution 
as an initial-value problem. For a single time step this explicit scheme scales like N 2 , with N denoting the 
number of grid points on a single time sheet. 


3.1. Hermiticity under the direct scheme 


Lemma 6. The direct scheme Eqs. (25) to (28) conserves Hermiticity of the density matrix. 


Proof. It is sufficient to show conservation for one time step under the assumption of Hermiticity of the 
initial density matrix placed on time sheets j Q — 1/2 and j Q . This is done in straight-forward way by showing 
that, within the scheme, complex conjugation is equivalent to a transposition of the density matrix on 
time-sheets j a + 1/2 and j a + 1. 
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Using M± = —M± and relations, such as 


{D’+pn^'iS?'- 1 ™)}* 


Jo'iJo 2 _ n 3 0 \3o 2 n Jo\3o 2 

P 21, ,• i l.,- ,■ P21 ■ l., ,■ P21.v , l „• . 

Jx'iJy 1 ' 2 iJxiJy Jx'iJy' 2 iJxiJy | ^ Jx' 2 ’Jj/'’« 

I - £■ 


n j 0 ;jo-2 
— P21 / 1 ■ • 

2 ->3 y '\3x->3y 3 x ~ ~*l i3y'\3&->3y 


A. 


y 


A, 


Jo 2 ’Jo _ jo 2 ’Jo jo 2 ’Jo _ Jo 2 ’Jo 

Pl2,- ,• ; I l Pl2 ■ • l dl2 •,-'+! , Pl2, , .,■/_! ,• , 

Jx ijy iJ x 'iJ y> \ 2 JxiJyiJx' iJy' 2 ^ JsiJyUjjT 2 Jy' JxiJy,J x 2 ’Jy' 


Ay 

= D'_ Pl2 (g[ jo - 1/2) -,g[ joY ), 

as well as the short-hand grid notation introduced above one may write 

1 


A t 


(29) 


/on(af o+i) ';af o+i) }* = 


J_ _ (M-h(^ Jo)/ )*—m + (g^ Jo) )*) La 

A 0 2 






j_ _ (-M + (gp o)/ )+Af + (g^ o) )) La 

A g 2 

- - D+p^g^-g^-^') 

= Pn(g[ jo+h) ;g[ io+hY ), 




(30) 


di2(^ io+i) ';# +1 V 


l 


A 


i Wi ^ A ° 

52F- : 

£>> 2 2 (^ o) '; ^ io) ) - D +Pll {g^ Y -g { 2 

l 


r — - 


i {-M + (g[ joY ) + M_{g ( ^ ) )) (jo -|) (jo) ' 

1_ _ (-M + (g< j ° ) ')+M_(g^° + ^ ) )) LA g 2 12 1 > 2 


A 0 


- £>-d22(# } ;^ o) ') - D_ Pll {g^ 0+h) -g^ +hY ) 


P2i(Gi j ° +1Y ,g i 1 jo+hY ) 


(31) 


/02i(^ Jo+1) ';# +l) )* = 


1 


1 M_(^° + ’ ) ) , -M + (^' e) ) 


i Af_(d^° + 2 ) )' -M+(g[ M f LA n 
A 0 2 


/02i(^ io) ';# i} ) 


- D'_pu(g?° +iy -,g[ jo+i) ) - D-p22(g¥ oY -,Q¥ o) ) 


i 


J_ _ -M_(g^° + ^')+M + (gt’°’) 
A 0 2 


1 M_(^ 0+|) )' + M + (^ ) ) 

(^LA- + 2 Pl2{Gl ’ G * } 


D' + p 11 (Q[ 3o ^ , ;Q ^ > ) - H+P22(ap o) ;^ o) ') 


■»(jo+ 2 ) . /->(jo + 2 
*1 » ^2 

Pi2(# +i) ;^° +1) ') , 


(32) 
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and 


P22(Q¥ 0+1Y \S¥° +1) ) 


1 

A a 


' 1 

-M_(g ( 2 Jo+ ^) La o 


M _ ( ^ + i)') _ M _ ( g^ + =)) UoY Uo) 

-X- P22{y 2 1*2 J 


^-Pi2(# +i) ';^' 0+1) ) - D + P2 i(^- 
1 


(io + 1 ) , ; gOo+J ) ) 


1 

A 0 


-M_(Sp + 2 ) ')+M_(Sp + = ) ) LA q 
2 

O’o+i). ^»0o+§) ^ pj ^ . /^Oo+i) 7 > 


1 Uo) _ (j 0 y 

-> -- P 22{y 2 1*2 ) 


D' +P21 (g^ J -,gr "') - d-pi 2 (02 ot *';S2 

= P22{Q^ +1) -A io+iy ) • 


(33) 


□ 


3.2. Stability of the direct scheme 

We first consider mass and potential to be position independent. Time-dependence is permitted. 
Proposition 7. For constant-in-position potential V and mass m, the direct scheme Def. [6] is stable for 


r 2 -p r 2 < _ 
x ' y — 4 


Proof. The proof is given by von Neumann analysis setting, for time-step A t , 


p+=e- A V k V k ' r 'p-, i,j = 1,2 


(34) 


(35) 


and inserting this expression into the scheme Def. [ 6 ] The resulting four linear homogeneous equations in 
of have a characteristic determinant which has the solutions 

' L J 


. 2 ,wA t .. P + p’ + ^~ . \/(P + P' + ^-) 2 - (1 + P 2 ){p ~ P 1 ) 2 
sin(-z-)|± = —- kk^± 


1 + / 


1 + /? 2 


(36) 


with 


;,( ) A £■( ) A 

= r 2 y sin 2 (^-A) + rl sin 2 (—— 


and 0 = Noting that 0,p,p' > 0, the maximum of the rhs of Eq. (36 1 is obtained for p = p' = 

Pmax = r 2 + r 2 and the positive sign for the square root. This leads to the condition 4 p m ax < 1 rendering 
sin 2 () | ± < 1. Finally, in the domain 0 < p^ < \. the radicant in Eq. (36) remains positive definite for 
0 > 0. □ 


Comments: 

(i) Von Neumann analysis reveals the single-cone energy-momentum dispersion relation of the scheme. 
For Hermitian p, one has k' = —k and p = p'. Hence for r 2 + r 2 < | the four real solutions for w are 
wi = oj 2 = 0 and 


lo 3 = —w 4 = — arcsm ■ 
A t 


/ 4p + 0 2 
1 + B 2 


(37) 


The origin of these four solutions can be interpreted as follows. For given k vector there is a positive- 
particle) and a negative-energy (antiparticle) solution. Both of them are needed for the construction of a 
complete basis set. Hence there are four types of matrix elements in Eq. (351: particle-particle (uq = 0) 
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Figure 1: Exact energy-momentum dispersion (meshed-surface) versus the direct-scheme dispersion: mass m = 0.02eV, 

v = 6.2 x 10 6 m/s, At = A = A x = A y , (a) A = 6nm (150 x 150 grid points), A t = 3.4 x 10~ 15 s (b) A = 0.6nm 

(1500 x 1500 grid points), A t = 3.4 x 10- 16 s. 

and antiparticle-antiparticle (012 = 0), particle-antiparticle (014) and antiparticle-particle (003). Thus the 
energy-momentum dispersion of the direct scheme is given by 


Ek = ± — arcsin 

A t 



(38) 


A comparison between the exact dispersion to the one for the direct scheme is given in Fig. |T] for a specific 
numerical example with m = 0.02eV. Note that for the case of 1500 grid points, the finite mass gap cannot 
be resolved for the scale used to allow the display of the entire range of k-vector space. 

(ii) Under the conditions of Prop. [7j preservation of positivity is ensured. Starting from a Hermitian 
positive density matrix p ~-, the latter may be expanded in the stationary orthonormal plane-wave eigenkets 


of the scheme according to Eq. (18) which evolve into according to Eq. (35) and k' = —k . 

We now turn to the case of time- and position- dependent mass and potential. The direct scheme cannot 
be expected to conserve positivity in general since it is based on a finite difference approximation of the 
linearized differential equation 


p+ = (l 


iA t H 

h 


)p~( 1 


i^tH 

h 


) = p 




2 

Hp~H 


P 


h 


[h, p~] 


Clearly, the original scheme is of Kraus form (with non-unitary Kraus operator), however, the linearized 
form is not. The maximum violation of positivity of the latter is of order (A t ) 2 ((^y) 2 for m± = 0) per time 

step. For states |$) for which (<b|p~|<l>) = 0 it is given by — (yy) 2 {&\Hp~H |<f>) (< 0 for positive p~). So 
starting from a positive definite initial density operator, it is conceivable that with time one looses positivity. 
Nevertheless, the direct scheme supports a conserved functional which replaces conservation of the trace for 
the continuum Dirac equation and which is closely related to the conserved functional of Lemma [3] for the 
bra-ket scheme. 
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Lemma 8. Under periodic, respectively, zero boundary conditions, 


Trip) = Tr {pn(^” 1/2) ; ^“” 1/2) ) + pn(G, 




A 0 D + p2i{g[ Jo) -,gi° 1/2) ) + A 0 D-p 12 (g.2° 1 / 2 ) ;g 2 o) )}, r = (r x = A 0 /A x ,r y = A 0 /A y \. 39) 

is a conserved quantity under the direct scheme. 

Here trace Tr, consistent with Definition |TJ implies the sum over (j x , j y ) £ Z 2 U (Z + |) 2 performed on 
the respective time-sheets]^] This conservation law corresponds to the conservation of the trace of p in the 
continuum limit, in which all spatial derivatives are carried out on a “single” time-sheet. 

Proof. Making use of the Cor. [l] (integration by parts on the lattice) one may write 


Tr 

E 

jx ijy 

E 

jx ijy 


{ P u(g? o+1/2) ;g[ jo+1/2) )+ P 22(g 2 


(jo + l) . n (jo + l) 


)} 


Pii(# ,+1/2) ;# +1/2) ) + P22(g : 


Ua+D.gUo + l^ 


Pii(g[ j °- 1/2) -,g[ j °- 1/2) ) - A 0 D +P21 {g^-g[ ia - 1/2) ) 


A 0 D'_ Pl 2 {g[ 3 o - 1 / 2 ) -g[ jo) ) 


E 

jx ijy 


P22(g { 2 jo) -,g { 2 jo) )- A 0 D_ Pl2 {g. 


G'o+i/2). ^(jo+i) 




A 0 D' + p 21 {g^ + 1 ) -gi 3o+1/2) ) 


= E 


Pii{g[ jo ~ 1 / 2 ) -,g[ jo ~ 1/2) ) - a 0 d +P21 (^ ) ;^- 1/2) ) 


A 0 D_pi2(5 2 ‘ 


Uo-l/2 ). r Uoh 


> »2 


E 

jx ijy 


P22(g 2 o) ',g 2 o) ) - A o D_p 12 (02 


O'o+l). ntio+l/2)} 


(Jo+l/2).g0'o+l)\ 


+ A 0 D + p 21 (t/)' 7 °~” , C/j 


= E 


Jx ijy 


pn(g^- 1/2) ;g^- 1/2) )+P22(g^W^) 


A 0 D +P 2 l {g ( f o) -g[ ] °- 1/2) ) + a 0 D_p 12 (g { 2 Jo ~ 1 / 2 ) -,g ( 2 Jo) ) 

A 0 D + p 21 {g[ Jo+ 1 ) -,g[ J ° +1/2) ) - A 0 D_p 12 (^ o+1/2) ;^ o+1) ) 


(40) 


□ 

In the time regime of positivity sufficient convergence conditions can be established from the conservation 


law Eq. (39). 


1 Note that p also is a matrix in spin space indicated by the subscripts (i,j) in p,,. 
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Proposition 9. For r x + r y < 1/2, with r„ = A a /A v , v £ {x, y}, and the assumption of positivity for p 
the scheme satisfies the estimate: 


Tr{pu + P 22 } < 


rpo 


1 - 2 (r x + r y ) 


(41) 


Proof. Since Hermiticity is already shown, we can write (here complex conjugation is denoted by the bar 
symbol~ for compactness of notation) 


_ l Y.iPiniUi J2iPi u i v i \ 

P V PiPiUi Y.iPiViVi ) 

for each time step j t £ (Z U Z + |). Under the assumption of positivity of p one has 0 < Pi < 1, real. We 
use the short-hand notation Tr{|M|} = j ) l-^1> Cor. IT1 as well as 2|5ft(it*u)| < ||w|| 2 + |M| 2 , to estimate 


followed by 


Tr {| - A 0 D+p 2 i + A 0 D_pi 2 |} 

= Tr{\-A 0 D + p 21 + A 0 {D' + p 21 }*\} 


= Tr{ | 

— A a D + p 2 i — A a {D + p 2 i} 1} 


= Tr { 25J (—A 0 D_|_p 2 i) } 




= 2 Tr < 

1 

Y.PiV V i + Ui 

. i 

-PiV V i Ui 

+ r x iR 

•OS 

_1 

< 2 Tr < 

r v \*l 

. i 

\+r v \M 

E k"? 

i 

Ui 

+ 2 Tr< 


T.Pi'F* 

_ i 

1 + r x \$s 

\ ^ X 

2^Pi y i 

. i 

-Ui 

< 2(r y + r x )T 

rhZviUiU 

+y ^PiVjVi > 


< 2(r y +r x )Tr{pn +P22}, 




T° = Tr {pn+P22 - A 0 L» + p2i + A 0 D_pi 2 } 

> Tr {pu + P 22 - | A 0 D+p 2 i - A 0 £)_pi 2 |} 

> Tr {pu + P22} — 2(r x + r y )Tr {pu + P22} ■ 


(42) 


(43) 


In Eq. (42) we denote a single-site shift of a spinor component by a superscript. For example, v/* denotes 
a shift of Vi by ±A y . □ 


In fact, a comparison of the conserved functional for the bra-ket and the direct scheme immediately 
allows a proof of Prop. [9] via Prop. [4j One observes that the difference between Tr{p} (Def. |T|) and is 
twice the difference between Tr{p} and E ° 


rpO 

r 


Tr{ Pll (#- 1/2) ;^- 1/2) ) + P 22(G i 2 jo) ;g { 2 jo) ) + 2A 0 K[U_ Pl2 (#" 1/2) ;# ) )]} (44) 

E° r + Ajk[D_p 12 {g^- 1/2) ] g^ ) )\ . 


This accounts for the factor of 2 in the denominator of the upper bound estimate in Eq. (41) as compared 
to Eq. (19). This observation immediately leads to an improved stability condition for the direct scheme 
via Prop. |5| 
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Proposition 10. Let r = r x = r y 

this scheme is stable and satisfies 


l/(2y/2) hold in (25 1 -(28 1 . Then, for all times of positive definite p 


Tr^l 1/2 ) + Tr (p&) < 2 T° , with T° = T° (2V5)i1/(2V2) • (45) 

Proof. Observe that T° (aV5)il/(2V5) = E 0 ^^ and use Pro P- 0 D 

Note that this estimate agrees with Prop. [7]for position-independent mass and potential term. Numerical 
simulations have confirmed this estimate, in particular, exceeding r = r x = r y = 1/(2V2) has led to 
instability in specific simulations. 

4. Numerical simulations 

(2+l)D Dirac fermions have recently been realized as low-energy excitations in graphene and topological 
insulators. In generic single-layer graphene there are four degenerate cones, while topological insulator states 
have an odd number of Dirac cones for a given surface. Although our simulations are rather generic, we use 
parameters typical for TIs, such as Bi 2 Te 3 , using a Fermi velocity v = 6.2 x 10 5 m/s. Simulation regions 
typically are lOOOxlOOOnm using at least 100x100 grid points per time sheet. At the rectangular boundary 
an absorbing layer is created by using a small imaginary contribution to the scalar potential which vanishes 
exponentially into the simulation region. The schemes were implemented in a MATLAB code using periodic 
boundary conditions allowing efficient global execution of the finite difference operations. 




Figure 2: Time-evolution of the single-particle density of an incoherent superposition of four Gaussian wave packets. 
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4-1. Gaussian superposition states 


We consider a mixed state of a zero-mass Dirac fermion consisting of four Gaussian wave packets centered, 
respectively, at = (±k x ,±k y ) with E x = h\k x \v = E y = h\k y \v = O.OleV. All initial Gaussians share the 
same spectral width, average energy, and center position, located at the center of the simulation region. For 
this simulation the bra-ket scheme is used to study the propagation under free-particle conditions (V = 0 
in the simulation region). Absorbing boundary conditions are implemented by means of an absorption layer 
via an imaginary potential contribution. m 

Figs. l(a)-l(d) show snapshots of the single-particle density of states at four different times. Compared to 
a Gaussian wave packet evolution under the Schrodinger equation, zero-mass Dirac fermions move dispersion- 
free. More precisely, due to the linear dispersion, plane-wave contributions move with Ev¥. While the 
present lattice schemes maintain a monotonic energy dispersion in form of a single cone on the lattice, the 
linear dispersion regime is limited to small values of k and A t , see Eq. (38) for the direct scheme. The 
lattice cone is represented by two energy bands of finite width which touch at k = 0 and display particle-hole 
symmetry. The analytic form for the dispersion underlying the bra-ket scheme was given in Ref. [221 . 

On a TI surface the number of flavors (Dirac cones) is odd and Dirac fermions are helical. Hence charge 
transport is linked to spin transport, with the particle spin being locked perpendicular to momentum. For 
the present simulation the spin current density is shown in several snapshots (corresponding to Fig. [2j of the 
expanding wave packets given in Fig. [3] For each snapshot the length of the ’’spin vector” is taken relative 
to the maximum value present at this particular time. For the superposition of four Gaussians, the initial 
state is spin-unpolarized by symmetry. However as the density current density evolves local spin density 
builds up and propagates along with the wave packets. These results have been obtained within the direct 
scheme as well. 








Figure 3: Time-evolution of the spin current density of an incoherent superposition of four Gaussian wave packets. 
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4-2. Particle injection from an electric contact 

For a second example we consider particle injection from an electric contact which is attached to one 
end of the rectangular simulation region. Next to an absorbing layer placed around the simulation region 
to absorb out-flowing contributions, particle injection from the contact is simulated using the linearity 
of the Dirac equation which is preserved by the difference schemes: we cut the density matrix into two 
orthogonal pieces according to position: the density matrix on the grid within the simulation region and 
the density matrix on the grid in the contact region. We first construct the density matrix from plane- 
wave solutions of the free-particle Dirac equation with occupation probabilities given by the Fermi-Dirac 
distribution function at temperature T = OK. We consider a filled Fermi sea of negative-energy states 
E < 0 and study electron injection between in the energy interval [0,-Ep]. Since the simulation is based on 
periodic boundary conditions it is somewhat advantageous to use injection based on a one-sided Fermi-Dirac 
distribution, i.e. to omit outgoing states from the start rather than damping them out in the absorption 
layer. In a first step we compute one time-step for the evolution of the truncated density matrix in the 
contact and record the contribution to the density matrix which arises on the grid of the simulation region. 
From then on this contribution is added step-by-step to the calculation of the time evolution of the density 
matrix on the grid of the simulation region. The basic assumption, common in transport simulations, is 
that particle injection provided by the contact is independent of the outflux of particles. 

The buildup of particle density due to the presence of a contact at the left end of the simulation region 
is shown in Fig. [4j For this simulation 120x120 mesh points, including 2 absorbing layers, adjacent and 
opposite to the contact as shown in Fig. (5^a), were used. A Fermi energy Ep = O.OleV and 50 randomly 
generated incident plane waves were used to construct the equilibrium density matrix. The simulation was 
performed over 500 time steps (third images in Fig. [4|. At final simulation time, steady-state is not quite 
reached. Also noticeable are two slight ridges in particle density propagating into the system. These can be 
attributed to the fact, that the construction of the equilibrium density matrix is based on plane waves for 
an infinite system, while the simulation is based on periodic boundary conditions (in y-direction). They are 
not resolved in the spin texture. A more detailed study of contact simulations based on this finite-difference 
approach will be given elsewhere. |25j 




Figure 4: (a) Charge buildup due to particle injection from a left-sided contact (particle reservoir) in thermal equilibrium at 
T = 0. The simulation region of 900nmx900nm is discretized by 120x120 mesh points. The corresponding applied bias is 
0.02eV. (b) Evolution of spin texture due to particle injection. 
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Examples for imaginary potential contributions for the modeling of absorbing boundary conditions are 
given in Fig. [5] 



Figure 5: (a) Imaginary potential contribution used to model absorbing boundary conditions: (a) absorbing layer at the contact 
end and opposite side only, (b) absorbing layer surrounding the simulation region. 






(d) 


Figure 6: Side view of the charge buildup due to particle injection from a left-sided contact (particle reservoir) in thermal 
equilibrium at T = 0. The simulation region of L x x Ly=900nmx900nm is discretized by 130x130 mesh points. The applied 
bias is 0.02eV. (a) 50 time steps, (b) 200 time steps, (c) 400 time steps, (d) 600 time steps. 


In Fig. [6]we show a side view of the charge build-up due to the action of an electric contact. Input data 
are as above, however, a grid size of 130x130 and the imaginary potential displayed in Fig. §b) was used. 
The time step A t is 3.9xl0 _15 s. Towards the central corridor of the simulation region fairly homogeneous 
charge distribution is reached. Due to the additional absorption regions (compared to Fig. [^a)) near y = 0 
and y = L y parallel to the average propagation direction x some loss of particle density is obtained at the 
y-boundaries. Steady state is reached at around 600 time steps. 

5. Extensions of the schemes 

5.1. External magnetic fields 

The numerical schemes presented here can directly be used to model the surface-state contributions to 
electron transport in topological insulators. These contributions inherently provide spin-polarized electric 
currents which may find application in spintronic devices. Disorder effects within this approach can be 
treated at a microscopic level via an explicit account of potential fluctuations. For this purpose it will be 
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useful to implement general electromagnetic potentials into the scheme. In presence of the vector potential 
the kinetic momentum entering the Hamiltonian Eq. ([2]) is modified from p to (p — ^A(x, y, <))rl and is 
incorporated via the Peierls substitution for the spinor components if on the lattice [21 EH [ 29 ]. For the 
present case it takes the form 


Pij 


U\tj 


Pij 


U;tj 

Xi,yi',x j,yj 


exp {-i(a% yz - ■ 


x j ,Vj 


)}p 


U;tj 

’■Jx i ,y i -x j ,y j 


(46) 


The real phase contribution off. is defined as the line integral of the vector potential A. starting at an 
arbitrary, but fixed position ( x 0 ,y 0 ) on the grid and ending on the lattice point ( Xi,yi ), 


q r( x ’V'> 

a xi,yi IT ds-A(s,t) | x=Xi,y—yi,t=ti ■ 

J ( x 0 ,y 0 ) 

Replacing p^ by the Peierls-transformed pij in Def. [4] and Def. [6] gives the respective numerical scheme 
in presence of an electromagnetic vector potential. Therefore, the gauge-invariant introduction of the elec¬ 
tromagnetic vector potential via the Peierls substitution on the grid leaves the stability properties of the 
respective scheme intact. 


5.2. Lindblad equation and Green’s function approach 

The two schemes presented here can be extended to a numerical treatment of the Lindblad equation for 
the Dirac Hamiltonian 

1 . ' (47) 


P(t) = ^ [H(t),p{t)] + D{p}, 


with the dissipator 




L„p{t)L+--{L+L^p{t)} 


(48) 


and initial condition p(0) = po , the dissipation term containing Lindblad operators = 1,... with the 

curly bracket {...,...} denoting the anti-commutator. Discretization of the dissipator follows the placement 
of the density matrix elements on the staggered grid, as discussed in Sect. |1.3| The Lindblad extension allows 
a treatment of Dirac fermions as an open quantum system, i.e. an account of dissipative effects (within a 
Markov approximation) in the system arising, for example, from electric contacts, phonons, nuclear magnetic 
moments, and particle transfer between surface and bulk states. 

Next to the density matrix approach, the Green’s function approach is the most versatile approach used 
in non-equilibrium electron transport theory. [5] Since the (single-particle) Green’s function fulfills the Dirac 
equation, in a real-space-time formulation of the latter a single-cone implementation on a lattice is ensured 
following the procedure outlined in this paper. The main difference is that the (single-particle) Green’s 
function inherently uses two time indices, while the density matrix only needs one (in the continuum limit). 
Hence, on the lattice introduced above, a single-particle Green’s function matrix element needs up to four 
time sheets. 


6. Summary and Outlook 

In summary we have proposed two closely related single-cone finite-difference schemes for the numerical 
treatment of the von Neumann equation for (2+l)D Dirac fermions. Main mathematical properties, such as 
stability the underlying energy-momentum dispersion relation, have been explored. Elementary numerical 
examples for Dirac fermion propagation characterized by particle density and spin texture have been given. 

Numerical solution of the time-dependent Dirac equation is quite memory intensive. Comparison of the 
two schemes shows that the advantage of the direct scheme regarding the number of matrix multiplications 


5 q is the fermion electric charge. 
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per time step is compensated by the reduced radius of convergence when compared to the bra-ket scheme. 
In MATLAB which we used for numerical analysis global matrix operations as compared to do-loops are 
most efficient. Therefore it is important to keep the number of full-sized matrices to a minimum. Our 
schemes are direct and are highly efficient since they scale proportional to the square of the number of 
spatial grid points. Nevertheless, the computation of the coherent part of the time-evolution is redundant 
since it is performed on left- and right-sided matrix indices independently. As evident from the bra-ket 
scheme, one may proceed as follows: decompose the initial density matrix in its orthonormal basis according 
to Eq. (18) and propagate the ket, transpose the ket to get the bra’s time-evolution, and construct the 
density matrix. [22] Depending on the number of basis functions needed to represent the problem this may 
be significantly less memory expensive. Extension to a treatment of the Lindblad equation along these 
lines, termed the quantum-jump approach, has been explored for the case of the Schrodinger equation. 126) 
Conceptual advantages of a density matrix or Green’s function formulation arise when dissipative effects and 
interaction processes, in general, need to be included, in particular for systems with a continuous spectrum 
typical to quantum transport. [30j 

Open questions arise regarding the treatment of simulation boundaries. In the present work a combi¬ 
nation of periodic and absorbing boundary conditions has been used to suppress communication between 
opposite simulation boundaries. We are currently exploring the potential of Lindblad operators to imprint 
the signature of electric contacts onto the simulation region. |25j Finally there are open questions regard¬ 
ing conservation of positivity of the direct scheme for the case of general position- and time-dependent 
electromagnetic potentials. 
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Appendix A. Explicit form of the direct scheme 


Here we give the explicit form of (25) - (28) for the direct scheme. Pairs of subscripts denote spatial 


grid position, the superscript denotes the time sheet. For clarity, we use double time superscripts throughout. 
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